Uncovering Key Metabolic Determinants of the Drug Interactions Between Trimethoprim and Erythromycin in Escherichia coli

Understanding interactions between antibiotics used in combination is an important theme in microbiology. Using the interactions between the antifolate drug trimethoprim and the ribosome-targeting antibiotic erythromycin in Escherichia coli as a model, we applied a transcriptomic approach for dissecting interactions between two antibiotics with different modes of action. When trimethoprim and erythromycin were combined, the transcriptional response of genes from the sulfate reduction pathway deviated from the dominant effect of trimethoprim on the transcriptome. We successfully altered the drug interaction from additivity to suppression by increasing the sulfate level in the growth environment and identified sulfate reduction as an important metabolic determinant that shapes the interaction between the two drugs. Our work highlights the potential of using prioritization of gene expression patterns as a tool for identifying key metabolic determinants that shape drug-drug interactions. We further demonstrated that the sigma factor-binding protein gene crl shapes the interactions between the two antibiotics, which provides a rare example of how naturally occurring variations between strains of the same bacterial species can sometimes generate very different drug interactions.


INTRODUCTION
A growing volume of microbiology research has been dedicated to characterizing drug-drug interactions in combination antibiotic treatment, because understanding how and why antibiotics interact with each other has the potential to help clinicians improve treatment efficacy, reduce dosage-related side effects and slow the emergence of antibiotic resistance (Tamma et al., 2012;Worthington and Melander, 2013;Mitosch and Bollenbach, 2014;Bollenbach, 2015;Wood, 2016). The use of high-throughput assays for measuring bacterial growth rate has accelerated the identification and characterization of antibiotic interactions (Yeh et al., 2006;Wood et al., 2012;Ocampo et al., 2014;Chevereau and Bollenbach, 2015;Zhou et al., 2015;Nguyen et al., 2016). Pairwise interactions between antibiotics can be classified as additive, synergistic or antagonistic ( Figure 1A) based on whether their combined effects on bacterial growth inhibition are equal to, greater than or less than those of their individual effects, respectively (Loewe, 1928(Loewe, , 1953Bliss, 1939;Yeh et al., 2006;Foucquier and Guedj, 2015). Suppression is a special type of antagonism in which the combined growth inhibitory effect is weaker than that of one antibiotic or both antibiotics alone (Yeh et al., 2006).
Our current understanding of the origins of antibiotic interactions still lacks comprehensive mechanistic underpinnings. There are several well-studied examples that explain why drugs that target the same pathways tend to interact synergistically. For example, trimethoprim and sulfamethoxazole inhibit two steps in the tetrahydrofolic acid synthesis pathway. Their synergistic interaction has been attributed to the drugs' ability to potentiate each other's inhibitory effects on key metabolite synthesis within this pathway (Minato et al., 2018). The synergism between two ribosome-targeting antibiotics of the macrolide family, lankacidin and lankamycin, can be explained by their ability to target neighboring sites in the ribosome (Auerbach et al., 2010). For antibiotics that target unrelated pathways, the ability of one drug to affect the uptake or efflux of another has been proposed to be a mechanism that results in drug interactions (Cokol et al., 2011;Brochado et al., 2018). For example, the synergism between streptomycin and penicillin was proposed to be due to the cell membrane-damaging effect of penicillin, which in turn facilitates bacterial cell penetration by streptomycin (Plotz and Davis, 1962). Conversely, antagonism can emerge when one drug limits the proton motive forcedependent uptake of another drug or increases its efflux through efflux pumps, resulting in decreased intracellular concentration of the latter (Brochado et al., 2018).
The suppressive effect of spiramycin, a ribosome-targeting antibiotic, on trimethoprim has been attributed to non-optimal regulation of ribosomal genes, which leads to an imbalance between cellular DNA and protein contents in the presence of trimethoprim-induced DNA stress (Bollenbach et al., 2009). This mechanism also underlies the antagonism between ciprofloxacin and tetracycline and has been proposed to hold more generally for antibiotic pairs comprising a DNA synthesis inhibitor and a protein synthesis inhibitor (Chait et al., 2007;Bollenbach et al., 2009). More broadly, polysaccharide and ATP synthesis in E. coli are key cellular processes that shape pairwise drug interactions between six antibiotics with different modes of action; supplementing growth medium with small molecule adjuvants that target these cellular functions alters drug-drug interactions in a predictable way (Chevereau and Bollenbach, 2015).
Erythromycin shares similarities with spiramycin with respect to their inhibitory effect on translation elongation (Kirst and Sides, 1989). Both are macrolide-type ribosome-targeting bacteriostatic antibiotics that bind to the 50S large ribosomal subunit, block nascent polypeptides (Hansen et al., 2002), and trigger the cold shock response in E. coli (van Bogelen and Neidhardt, 1990). Erythromycin is expected to suppress trimethoprim. However, the interaction between erythromycin and trimethoprim was experimentally characterized as additive (Ocampo et al., 2014) or predicted to be additive against E. coli (Wood et al., 2012). What causes this apparent lack of a suppressive interaction between erythromycin and trimethoprim? To address this question, we applied an RNA-Seq transcriptomics approach to capture the effects of single and combined antibiotic treatments on gene expression in E. coli. Gene regulation responses often show conflicts when antibiotics that generate different transcriptional responses are combined. Bacteria often resolve such conflicts by prioritizing gene expression to one of the two drugs or by averaging the gene expression due to both drugs (Bollenbach and Kishony, 2011;Young and Elowitz, 2011). We explored whether this phenomenon could help us identify genetic determinants that shape the interactions between erythromycin and trimethoprim. Notably, our study revealed that when the two antibiotics are combined, the expression of sulfate reduction genes escaped the stronger overall effect of trimethoprim on the transcriptome. This unexpected transcriptomic response underlies the lack of a suppressive interaction in a sulfate-limited growth environment. In addition, we showed that the presence of an intact sigma factor-binding protein gene crl also contributes to a suppressive interaction between the two antibiotics. Overall, our findings highlight the roles of gene expression prioritization patterns of an essential metabolic pathway and a key genetic background variation in shaping drug interactions.

Growth Medium, Antibiotics and Bacterial Strains
LB broth Miller (pH 7.4) was used as the default growth medium throughout this study, and all chemicals were obtained from Sigma-Aldrich (St. Louis, United States). Where required, LB broth was supplemented with 2 mM sodium sulfate. Trimethoprim (TMP) was dissolved in dimethyl sulfoxide, while erythromycin (ERY) and spiramycin (SPR) were dissolved in ethanol to obtain 50 mg/mL stock solutions. The E. coli MG1655 strain in this study is a crl-deficient strain used in a previous work with a 777 bp IS1 transposable element insertion in the crl coding region (Bollenbach et al., 2009). The E. coli BW25113 ancestral strain and the crl single gene knockout derivative strain of BW25113 are from the Keio collection (Baba et al., 2006). To account for the kanamycin resistance marker that is present in the crl knockout strain, BW25113 was transformed with the low copy number pUA66 plasmid (Zaslaver et al., 2006) to generate the kanamycin-resistant BW25113/pUA66 control strain. The IS1 insertion in crl of our MG1655 strain and the intact crl in the BW25113 parental strain were verified using the primer pair 5 -CAGGAAATCACCGACTGGAT-3 and 5 -CGACGTCGGTGCTACGTATT-3 (Freddolino et al., 2012).

Growth Rate Measurement and Data Analysis
Prior to each kinetic growth rate experiment in 96-well microtiter plates, E. coli strains were pre-cultured overnight for approximately 20 h in antibiotic-free LB broth with shaking (225 rpm) at 37 • C. On the next day, antibiotic stock solutions for trimethoprim and one of the two macrolides (ERY or SPR) were diluted in LB broth to obtain linear FIGURE 1 | In contrast to another ribosome-targeting antibiotic spiramycin (SPR), erythromycin (ERY) does not suppress trimethoprim (TMP) when tested against a crl-deficient MG1655 strain. (A) Schematic of different drug-drug interaction isoboles in a two-dimensional drug concentration space. Loewe additivity of drug-drug interactions is used for quantifying effects of drug combinations. Lines of constant inhibition (isoboles) in the 2D drug concentration space represent combinations of drugs A and B that result in same normalized growth rate. The four types of interactions shown are synergistic (green), additive (red), antagonistic (orange) and suppressive (blue). (B) TMP and ERY show almost additive interaction against crl-deficient E. coli MG1655. Isoboles in grayscale represent exponential growth rates normalized to that of the no-antibiotic treatment control in 2D concentration space. When tested against a crl-deficient E. coli MG1655 strain, ERY and TMP show an additive interaction. In contrast, SPR exerted a mildly suppressive effect on TMP under the same experimental conditions. concentration gradients, which were subsequently used to set up two-dimensional antibiotic concentration gradients in clear, flat-bottom 96-well microtiter plates with lids (Thermo Fisher Scientific, United States). The overnight pre-culture was diluted 1:1,000 in 200 µl cultures with and without antibiotics. Six technical replicates of cultures with no antibiotic treatment were included in each microtiter plate as control. The microtiter plates were incubated for 16-18 h at 37 • C with continuous shaking on a Synergy H1 Hybrid Multi-Mode Reader (BioTek, United States). Absorbance at the 600 nm wavelength (OD 600 ) was measured every 10 min. Only growth curves of sub-lethal antibiotic treatments were retained for further analyses. Growth curves that showed lethal effects associated with high concentrations of single or combination antibiotic treatments were excluded even if an initial period of transient growth can be recorded.
Calculation of exponential growth rates was performed in MATLAB (MathWorks, United States). First, OD 600 measurements from wells that contained blank LB broth were subtracted from all data points. Blank-corrected absorbance measurements (OD 600 ) that lie within the early exponential phase (defined here as 0.02 < OD 600 < 0.18) were logtransformed. An array of linear regression values within the interval was calculated for each culture using all combinations of six consecutive data points for every bacterial culture. A 60-min sliding window determined the maximum exponential growth rate within this interval. Normalized exponential growth rates of bacterial cultures with antibiotic treatments were calculated by normalizing their maximum exponential growth rate to that of the no-antibiotic treatment control in the same experiment. Each growth rate experiment was performed three times to obtain three independent biological replicates. Arrays of average normalized exponential growth rates in the [0, 1] interval were thus obtained for each experimental condition involving different E. coli strains, antibiotic pairs or growth media, which were noise-filtered using the ordfilter2 function with a round filter. Isoboles representing the same levels of normalized exponential growth rates in the two-dimensional concentration space of antibiotic pairs (isocontours) were visualized using the contour function.

RNA-Seq
Prior to RNA isolation, bacterial cultures were grown in the same conditions as described above in LB medium (Miller). Two biological replicates of total RNA samples were isolated for each antibiotic-treatment and the no-treatment control ( Table 1). Eight technical replicates were included for every condition for each biological replicate. After 3 h ± 20 min of growth, exponential phase bacterial cultures (OD 600 < 0.35) were mixed with RNAprotect Bacteria Reagent (Qiagen, Germany), and total RNA was extracted using the SV Total RNA Isolation System (Promega, United States) according to the manufacturers' instructions with minor modifications. The RNA integrity numbers (RIN) of all the total RNA samples were assessed using the Agilent RNA 6000 Nano Kit and Agilent 2100 Bioanalyzer (Agilent Technology, United States) and were in the 8.8 ≤ RIN ≤ 10 range.
Illumina high-throughput sequencing was performed by the Vienna Biocenter Core Facility. Ribodepletion of total RNA was carried out using the Ribo-Zero rRNA Removal Kit (Epicenter, United States). mRNA library was generated using the NEBNext Ultra Library Prep Kit for Illumina (NEB, United States). For multiplexing, TruSeq adaptors were ligated to individual samples. Paired-end sequencing was carried on HiSeq 2500 V4 with read length of 125 bp (Illumina, United States). Combination treatment with trimethoprim and erythromycin that resulted in approximately 50% reduction in exponential growth rate relative to the no-antibiotic control group

Transcriptomic Data Analysis
Raw Illumina reads were trimmed and filtered using the NGS QC Toolkit (Patel and Jain, 2012). The 3 and 5 ends of raw reads were trimmed when Phred score < 20. Trimmed reads that were shorter than 125 bp were discarded. Next, filtered reads with Phred score > 20 throughout 80% of their lengths were retained, while reads with ambiguity in > 2% of all bases were discarded. Illumina adaptors were removed using TrimGalore (Martin, 2011). Using Bowtie 2, the trimmed and filtered reads were mapped to the E. coli MG1655 reference genome with gene annotations (genome assembly ASM584v2) from Ensembl Bacteria (Langmead and Salzberg, 2012;Yates et al., 2016). BAM/SAM file manipulations were carried out using SAMtools commands for viewing, sorting and indexing (Li et al., 2009). All BAM ( * .bam), BAM index ( * .bam.bai) and reference genome files are available on the Dryad repository (see Supplementary  Table 1). Gene counts were extracted from the sorted and indexed BAM files using the Python script HTSeq (Anders et al., 2014).
After excluding genes with no reads, the median sequencing coverage was 683x.
Using the DESeq2 Bioconductor package (Anders and Huber, 2010;Anders et al., 2013), genes with an average of fewer than 10 reads per sample across all the experimental conditions were discarded. For each gene, the geometric mean for gene counts was calculated across all sample groups, which was applied as the normalization factor for that gene. Within each sample, the median of all normalized gene counts was applied as the size factor for that sample in the second normalization step. For principal component analysis (PCA), logarithmic transformation of the normalized count data was carried out to remove the dependence of the variance on the mean. PCA was performed using default parameters of the plotPCA function. The PCA score plot was generated using the 500 most variable genes across all the sample groups in our dataset. Briefly, the covariance matrix computed from the logarithm of normalized count data was decomposed into eigenvectors, which were used as weighting coefficients to calculate the principal component scores (Todorov et al., 2018). Top/bottom gene loading on the PCA were calculated using the Bioconductor package pcaExplorer (Marini and Binder, 2019). Heat-map diagrams and dendrograms from the two-way hierarchical clustering analysis of the 100 most variable genes were generated using the heatmap.2 function. For every gene, the logarithm of normalized gene counts was transformed into Z-scores such that their mean is 0 and standard deviation is 1 across all the sample groups.
For differential gene expression analysis, a negative binomial general linear model was fitted for each gene, and the Wald test was applied for significance testing. In samples that were treated with one or two antibiotics, genes are classified as differentially expressed relative to the no-treatment control group if their adjusted P-values are less than 0.01 after controlling for false discovery rate using Benjamini-Hochberg correction (Benjamini and Hochberg, 1995). Differentially expressed genes were compiled for five pairwise comparisons: (1) TMP against nodrug control, (2) ERY against no-drug control, (3) TMP + ERY against no-drug control, (4) TMP + ERY against ERY and (5) TMP + ERY against TMP. Venn diagrams were constructed for the differentially expressed genes in pairwise comparisons (1), (2) and (3). Pathway enrichment analysis was performed for pairwise comparison (5) using the GAGE (Generally Applicable Gene-set Enrichment) package (Luo et al., 2009). Enrichment of the assimilatory sulfate reduction KEGG pathway (eco00920) in a wider context of sulfur metabolism was visualized using the Bioconductor package Pathview (Luo and Brouwer, 2013;Kanehisa et al., 2016). Across pairwise comparisons (1), (2) and (3), two-way hierarchical clustering analysis of average log 2 fold-change for cys genes was performed using the heatmap.2 function.

Quantitative Real-Time PCR
In three biological replicates, E. coli MG1655 was incubated at 37 • C overnight with shaking in 2 mL antibiotic-free LB medium as three independent cultures. On the next day, saturated precultures were diluted 1:500 in 2 mL fresh LB medium (Miller) containing TMP, ERY + TMP or ERY + TMP + 2 mM sodium sulfate, where ERY and TMP concentrations are 42 and 0.0625 µg/mL, respectively. After approximately 3 h of further incubation with shaking (225 rpm), total RNA was extracted from each bacterial pellet using a phenol-chloroformed based method with QIAzol Lysis Reagent (Qiagen, Germany). After treatment with DNase I solution at 37 • C (New England Biolabs), 200 ng total RNA of each sample was added as template for cDNA synthesis using the LunaScript RT SuperMix Kit (New England Biolabs, United States) according to the manufacturer's instructions.
For qRT-PCR, equal volumes of each cDNA template were amplified using Luna Universal qPCR Master Mix (New England Biolabs, United States) and 0.25 µM gene-specific oligonucleotide primer pairs for the target gene cysH (5 -CTGAATGCCAAAGCTGGAAG-3 ; 5 -AACGCCGAACTGGA AAAACT-3 ), as well as the internal reference genes rpoB (5 -GTAAGGCACAGTTCGGTGGT-3 ; 5 -ATTTCCTGCAGGG TGTATGC-3 ) and secA (5 -GGTAGTCGTAACGATCGCA-3 ; 5 -TTTCCATCTCCGGTTCCAT-3 ). The thermocycler profile on the LightCycle 480 System (Roche, United States) comprised 60 s of initial denaturation at 95 • C, followed by 45 cycles of denaturation (15 s at 95 • C), extension (30 s at 60 • C) and fluorescence signal acquisition. Reactions were carried out in two technical replicates for the three biological replicates. Using the 2 − CT method, the relative expression levels of cysH were normalized to the average expression of rpoB and secA. Batch effect was accounted for by normalizing cysH expression levels in the TMP + ERY and TMP + ERY + 2 mM sulfate groups to those in the TMP-only reference group within each set of biological replicates. Relative expression levels of cysH are expressed in log 2 fold-change. Error bars denote standard error of the mean.

RESULTS
To characterize the interactions between trimethoprim (TMP) and erythromycin (ERY) against an E. coli MG1655 wild-type strain, we set up a two-dimensional concentration gradient of the two antibiotics in LB medium and measured the optical density (600 nm) of bacterial cultures over time using kinetic assays. All experiments were carried out in LB medium at 37 • C due to the availability of existing literature that focused on the inhibitory effects of sub-lethal TMP or ERY treatments under similar experimental conditions (Chevereau and Bollenbach, 2015;Giroux et al., 2017;Chen and Wang, 2021). The maximum growth rates of bacteria were quantified during the early exponential phase. The isoboles in the contour plots represent pairwise combinations of antibiotic concentrations that resulted in the same exponential growth rates relative to the noantibiotic treatment controls. The overall shape of the contours suggests that the interaction between TMP and ERY is almost additive (Figure 1B). We then performed the same experiment for TMP and spiramycin (SPR) and found that SPR has a clearly antagonistic, and even mildly suppressive effect on TMP (Figure 1B), which is in agreement with previous results despite the differences in experimental conditions (Bollenbach et al., 2009). Our results confirmed that ERY and SPR show different pairwise interactions with TMP under the growth conditions in this study.
To identify genes and cellular pathways that shape the interaction between TMP and ERY, we obtained the RNA-Seq profiles of exponentially growing cultures treated with single antibiotics at their half-maximal inhibitory concentrations (IC 50 ), as well as with a combination of these TMP and ERY concentrations using paired-end Illumina HiSeq (Table 1). Principal component analysis (PCA) was used to dissect how each treatment contributes to the total variance in gene expression and visualize the distribution of sample groups on a PCA score plot (Figure 2A). TMP and ERY triggered distinct transcriptomic responses and formed separate clusters. The position of the TMP + ERY combination treatment group is highly similar to that of TMP on PC1, which explains most of the variance (78%). It is mostly on PC2, which captures 16% of the variance, where the TMP + ERY group differs from the TMP group. Differential gene expression analysis relative to the no-antibiotic control group was performed (Supplementary File 1). The Venn diagram in Supplementary Figure 1A shows that the number of differentially expressed genes (Wald test: P < 0.01 after Benjamini-Hochberg correction) was much higher for TMP (2,443 genes) compared to ERY (1,008 genes). Linear regression analysis (Supplementary Figure 1B) revealed a high degree of similarity between the log 2 fold-changes of differentially expressed genes that are common to the TMP and TMP + ERY groups (coefficient of determination, R 2 = 0.925). Overall, the changes in transcriptomic profiles of the TMP + ERY group are driven by the stronger transcriptomic responses associated with TMP.
Next, we performed a two-way hierarchical clustering analysis to contrast the gene expression profiles of the top 100 genes with the greatest variations across all sample groups. Almost all the genes in the TMP + ERY group displayed gene expression responses that are consistent with those in the TMP group ( Figure 2B). The only clear exceptions are the cysH gene and cysIJ operon from the assimilatory sulfate reduction pathway, which belongs to the sulfur metabolism gene ontology (GO) category. FIGURE 2 | In TMP + ERY combination treatment, TMP exerts a much stronger overall effect on the transcriptomic response than ERY does. (A) TMP triggers a different transcriptomic response than ERY but is similar to the combination treatment of TMP + ERY. Principal component 1 (PC1) explains 78% of the total variance in gene expression, while PC2 accounts for 16%. TMP (blue) and ERY (yellow) trigger different transcriptomic responses, and the response to TMP + ERY treatment (green) is more similar to TMP than it is to ERY. (B) Except for three sulfur metabolism related genes, TMP and TMP + ERY show a high degree of similarity in gene expression profiles. Normalized gene counts from RNA-Seq were Z-score transformed across all sample groups before a two-way hierarchical clustering analysis was performed based on the top 100 genes with the highest variance. The heatmap revealed that the transcriptomic responses for the TMP and TMP + ERY groups were highly similar across four MultiFun gene ontology (GO) categories for 97 genes. The only exception is the sulfur metabolism category in which cysH and the cysIJ operon do not display TMP-induced upregulation in the TMP + ERY group (inset).
The expression levels of these three genes in the TMP + ERY group resemble those in the ERY and the no-drug control groups. Differential gene expression analysis indicates that TMP led to significant upregulation of cysH and cysIJ relative to the noantibiotic control group, while the TMP + ERY and ERY groups showed expression levels that are similar to the control group (Figure 2B inset). In other words, ERY in the TMP + ERY combination treatment lowered the expression of cysH and cysIJ relative to that in the TMP-treatment group.
We then examined the wider implications of the distinct expression patterns for cysH and cysIJ. Sulfate reduction genes are known to be differentially expressed when E. coli is treated with TMP (Sangurdekar et al., 2011;Mitosch et al., 2017). The assimilatory sulfate reduction pathway ( Figure 3A) scavenges extracellular sulfate and reduces intracellular sulfate to sulfide via several intermediate steps (Sekowska et al., 2000). This pathway is activated by the transcriptional activator CysB and the N-acetyl-L-serine inducer (Ostrowski and Kredich, 1989; FIGURE 3 | Transcriptomic response of the assimilatory sulfate reduction pathway. (A) Genes and operons of the assimilatory sulfate reduction pathway. The underlined genes and operons were included in the two-way hierarchical clustering analysis. cysG was excluded from this analysis due to its different regulatory mechanism from all other genes in the sulfate reduction pathway. (B) The expression profiles of the sulfate reduction pathway follow the same trend as cysH and cysIJ. Two-way hierarchical clustering analysis for genes of the sulfate reduction pathway and their transcriptional regulator gene cysB. With the exceptions of cysB and sbp, all the cys genes in this pathway generally follow the same gene expression profiles as those of cysH and cysIJ across all the sample groups. (C) Average log 2 fold-change of sulfate reduction gene expression levels in the TMP, ERY and TMP + ERY groups relative to the no-drug control group. Neidhardt, 1996). The end-product of the pathway, sulfide, depletes the N-acetyl-L-serine inducer via the O-acetyl-L-serine precursor, and thus exerts negative feedback on the expression of cys genes (Neidhardt, 1996). CysB is also negatively regulated by itself and N-acetyl-L-serine (Jagura- Burdzy and Hulanicka, 1981;Neidhardt, 1996). We hypothesized that the distinct expression patterns of cysH and cysIJ in the TMP + ERY group extends to other sulfate reduction genes, because they share the same transcriptional regulatory mechanism.
We performed a separate two-way hierarchical clustering analysis for cysB and the sulfate reduction genes of Figure 3A using row Z-scores and average log 2 fold-change relative to the no-drug control. With the exception of cysB and the periplasmic sulfate-binding protein gene sbp, the expression profile of all the cys genes in the sulfate reduction pathway followed the same trend as cysH and cysIJ (Figures 3B,C). To exclude the possibility that this was a consequence of differences in growth rates across the four sample groups, we compared the TMP (IC 50 ) and ERY (IC 50 ) groups with an additional TMP + ERY sample group with approximately 50% reduction in exponential growth rate. This was achieved by reducing TMP and ERY concentrations for this combination treatment group. Two-way hierarchical clustering analysis confirmed that the TMP (0.2 × IC 50 ) + ERY (0.6 × IC 50 ) group clustered with ERY rather than with TMP ( Supplementary  Figure 3), which is consistent with the trend in Figure 3B. This suggests that the sulfate reduction genes are generally restrained from responding to TMP when ERY is also present.
Treating E. coli with TMP in LB medium causes oxidative stress (Giroux et al., 2017), which is expected to compromise sulfate reduction and promote auto-oxidation of sulfide (Dolla et al., 2006). LB medium is also considered to be a sulfate-poor growth medium and contains approximately 0.1-0.15 mM of sulfate (Kertesz, 2000). The upregulation of the sulfate reduction pathway in the presence of TMP is a transcriptomic signature of sulfur limitation and decrease in intracellular sulfide availability (Neidhardt, 1996;Gyaneshwar et al., 2005). Thus, it is consistent with the negative feedback transcription regulation for this pathway that counteracts sulfur limitation (Neidhardt, 1996). Further links between oxidative stress and upregulation of the sulfate reduction pathway had previously been established by transcriptomic studies in which E. coli cultures grown in LB medium were treated with oxidizing agents such as chlorine, hydrogen peroxide, paraquat and sodium salicylate (Pomposiello et al., 2001;Wang et al., 2009). These are consistent with the TMP-induced upregulation of assimilatory sulfate reduction genes we observed. However, it is conceivable that other TMPinduced effects on one-carbon metabolism, as well as methionine and S-Adenosyl methionine (SAM) synthesis could also have contributed to the observed upregulation of the assimilatory sulfate reduction pathway.
Using the TMP versus no-drug control and TMP + ERY versus TMP pairwise comparisons, we performed gene set analysis on the assimilatory sulfate reduction pathway and neighboring nodes of the sulfur metabolism network (KEGG annotation). ERY in the combination treatment suppressed the upregulation of cys sulfate reduction genes when TMP + ERY was compared against TMP (Supplementary Figure 4B). A volcano plot for the TMP + ERY versus TMP comparison also confirmed that cys genes were amongst some of the most downregulated genes (Supplementary Figure 5). An intriguing observation that emerged from the pathway analysis was the upregulation of tauD and ssuE/ssuD, whose gene products allow taurine and alkanesulfonates, respectively, to be utilized as alternative sulfur sources under sulfate starvation (van der Ploeg et al., 1996;Eichhorn et al., 1999). The lack of an expected upregulation of the sulfate reduction pathway in the TMP + ERY group implies a transcriptional response to sulfur limitation that results in lower fitness.
To validate our key observation that ERY in the TMP + ERY combination treatment suppresses the expression of the sulfate reduction pathway, we performed quantitative real-time PCR (qRT-PCR) on E. coli MG1655 that had been treated with TMP and TMP + ERY. We focused on the PAPS reductase gene cysH, which showed one of the highest variations in expression levels across sample groups ( Figure 3B). cysH expression was normalized to the average expression of two stably expressed essential genes (according to our RNA-Seq based differential gene expression analysis), rpoB (encoding the RNA polymerase β-subunit) and secA (encoding the protein translocase subunit SecA) as internal reference genes. As expected, combining ERY with TMP resulted in an approximately twofold reduction in cysH expression level relative to the TMP-treated reference group ( Figure 4A). Under sulfur limitation, supplementing growth media with sulfate is expected to stimulate the expression of the sulfate reduction pathway Wheldrake and Pasternak, 1965). To verify this, we added 2 mM sodium sulfate to the LB medium that contained TMP + ERY treated MG1655 and observed an increase in average cysH relative expression compared to that without exogenous sulfate. This increase, however, did not fully offset the suppressing effect of ERY on cysH expression.
We hypothesized that a transcriptional response to sulfur limitation that is detrimental to fitness in the TMP + ERY group shapes the drug interaction between TMP and ERY. To test this hypothesis, we supplemented LB growth medium with 2 mM sodium sulfate and performed the same growth rate assays as before in TMP + ERY two-dimensional concentration space. Notably, ERY exerts a suppressive effect on TMP in the presence of additional exogenous sulfate, especially at higher TMP concentrations (Figure 4B). The disproportionately large beneficial effect of sulfate on normalized growth rate in the presence of both drugs altered the shape of the isoboles in the 2D drug concentration space (Supplementary Figure 6B).
To expand our characterization of TMP + ERY interactions to another key reference strain of E. coli, we performed additional growth rate assays using the BW25113 strain (Grenier et al., 2014). Although both BW25113 and MG1655 were derived from the E. coli K12 ancestral strain, we observed that ERY exerted a pronounced suppressive effect on TMP in the two-dimensional contour plot for BW25113 ( Figure 5A). A major difference between BW25113 and the crl-deficient MG1655 strain used in this study is that the former carries an intact crl gene. Crl is a sigma factor-binding protein and transcription activator that exerts global regulatory effects on the transcriptome by increasing the affinity of RNA polymerase for the principal sigma factor of the stationary phase σ S (Typas et al., 2007;Cavaliere and Norel, 2016;Cartagena et al., 2019;Xu et al., 2019). The 86 genes that are controlled by Crl in E. coli K-12 form 40% of the RpoS regulon (Santos-Zavaleta et al., 2019). The loss of crl can be detected in E. coli strains under laboratory settings that favored the inactivation of the general stress response (Faure et al., 2004), as well as in naturally-occurring pathogenic strains of E. coli (Provence and Curtiss, 1992;Asadi et al., 2018). We repeated the growth rate experiment using the BW25113 crl knockout mutant from the Keio collection to mimic the effects of crl-inactivation in the MG1655 strain we used. As expected, deleting crl abolished the suppressive drug interaction previously observed (Figure 5B), yielding an antagonistic interaction. Deleting crl substantially altered the susceptibility of BW25113 to TMP. The IC 50 for TMP approximately doubled from 0.10 µg/mL to 0.21 µg/mL, while IC 50 for ERY only decreased marginally from 48 µg/mL to 42 µg/mL. On the vertical axes of the contour plots, it is evident that the change in shapes of isoboles was mainly driven by the increase in tolerance to TMP in the absence of crl (Figure 5). The remaining difference between TMP + ERY interactions in the MG1655 (additive; Figure 1B) and BW25113 crl strains (antagonistic; Figure 5B) is most likely to be due to other variations in their genomes.

DISCUSSION
Applying a transcriptomic approach to identify genetic determinants that affect the interaction between two antibiotics with distinct modes of action is challenging, because the transcriptomic signatures of individual antibiotic treatments are superimposed. We used the interactions between TMP and ERY as a model system, because ERY and SPR show different FIGURE 4 | Perturbing the sulfate reduction pathway by adding exogenous sulfate increases cysH expression level in TMP + ERY treated MG1655 and changes the TMP + ERY interaction from additive to suppressive. (A) qRT-PCR validation of cysH expression. Relative to the TMP-treated MG1655 reference group, the TMP + ERY treated group showed an approximately twofold reduction in cysH expression, which confirmed the role of ERY in suppressing cysH expression. Supplementing the LB growth medium in which the TMP + ERY treated group is grown with 2 mM sulfate caused an increase in average cysH relative expression. (B) ERY exerts suppression on TMP in growth medium supplemented with sulfate. Isoboles for TMP + ERY drug combinations when tested against the E. coli MG1655 strain. Supplementing LB medium with sodium sulfate at 2 mM changed the interaction between TMP and ERY from additive to suppressive. Suppression by ERY on TMP is stronger at higher concentrations of TMP in the presence of exogenous sulfate.
interactions with TMP ( Figure 1B) despite their similarities as ribosome-targeting antibiotics. TMP inhibits folic acid metabolism, depletes key metabolites and induces DNA damage (Sangurdekar et al., 2011;Giroux et al., 2017). Activating global stress responses such as the SOS response and the ppGppmediated stringent response has wide-ranging effects on the transcriptome (Courcelle et al., 2001;Khil and Camerini-Otero, 2002;Durfee et al., 2008;Traxler et al., 2008). Despite this, the transcriptomic response of the sulfate reduction pathway in the TMP + ERY group prioritized to that of ERY rather than to TMP (Figure 3 and Supplementary Figure 3). By supplementing the culture medium with exogenous sulfate and perturbing the affected pathway (Figure 4A), we successfully changed the additive interaction between TMP and ERY to a suppressive interaction ( Figure 4B). These results demonstrate that analyzing prioritization patterns of conflicting transcriptomic responses can be a useful tool for dissecting drug interaction mechanisms.
Given the disproportionately large beneficial effect of sulfate addition on relative exponential growth rates when the two drugs were combined, we argue that the suppressing effect of cys upregulation by ERY is a transcriptional response that does not maximize fitness. The sulfate reduction pathway is a subset of sulfur metabolism, which is subject to strong negative feedback control, because L-cysteine is one of the most toxic amino acids at high concentrations (Datta, 1967;Avcilar-Kucukgoze et al., 2016). In addition to the ability of sulfide to deplete O-acetyl-L-serine (precursor to the N-acetyl-L-serine inducer), negative feedback is also mediated by the inhibitory effect of L-cysteine on the conversion of serine to O-acetyl-L-serine (Neidhardt, 1996). TMP upregulated cys genes of the sulfate reduction pathway and slightly downregulated cysB (Figures 3B,C), and these changes are consistent with the transcriptional regulatory mechanism of this pathway in response to sulfur limitation. The sulfate binding protein gene sbp is the only sulfate reduction gene that showed FIGURE 5 | The suppressive effect of ERY on TMP requires the presence of the crl gene in BW25113. Isoboles for TMP + ERY drug combinations for the E. coli BW25113/pUA66 and BW25113 crl strains. Whereas ERY exerts a strong suppressive effect on TMP when tested against the BW25113/pUA66 strain which has an intact crl gene (A), the suppression is removed when tested against the BW25113 crl strain (B). similar expression levels in both the TMP and TMP + ERY groups ( Figure 4B), but sbp is known to be activated more specifically by sulfate starvation (van der Ploeg et al., 1996). Given how tightly sulfur metabolism is regulated, it is surprising that the cys genes displayed an unusual transcriptional response in TMP + ERY (Figures 3B, 4A). As a ribosome-targeting antibiotic, ERY is not expected to directly affect the intra-and extracellular availability of metabolites for sulfur assimilation. More metabolomic insights into the regulation of the cys genes in response to combination treatment may reveal why TMP + ERY treated cells displayed a transcriptional response that is unlikely to be effective in counteracting the effects of TMP in terms of oxidative stress and sulfate limitation.
We also demonstrated that the sigma factor-binding protein gene crl plays a significant role in determining the interaction between TMP and ERY against E. coli. This finding lends support to the growing recognition that interactions between the same combinations of antibiotics can sometimes differ between bacterial species, as well as between different strains of the same species (Brochado et al., 2018). For example, the interaction between nalidixic acid and tetracycline against several multidrug resistant clones of Acinetobacter baumannii and E. coli is synergistic, while their interactions in susceptible clones of the same species are not (Gaurav et al., 2021). A previous study showed that drug interactions between antibiotics representing six main modes of action are surprisingly robust against most single gene deletions in E. coli, and that only a very small minority of single gene knockout strains produced qualitative changes in drug interactions (Chevereau and Bollenbach, 2015). Therefore, crl represents a rare example of single gene deletions that caused such shifts in drug interactions. While it is not clear if there are any interactions between the Crl-RpoS regulon and the assimilatory sulfate reduction pathway, our findings hint at the possibility that the effects of specific gene expression patterns on drug interactions may be contingent on the presence or absence of other genes such as transcriptional regulators. Additional future efforts could also be directed toward systematic discovery of specific genetic contexts in which alterations in drug interactions arise due to variations between strains, including the gain of antimicrobial resistance genes and loss of transcriptional regulators.
In conclusion, drug interaction is a complex phenotype that is determined by antibiotic modes of action (Cokol et al., 2011;Ocampo et al., 2014;Chevereau and Bollenbach, 2015), gene expression patterns (Bollenbach et al., 2009;Bollenbach and Kishony, 2011), genetic variations across species and strains (Chevereau and Bollenbach, 2015;Brochado et al., 2018;Gaurav et al., 2021). Our study exemplifies that metabolite abundance in the growth environment and gene prioritization patterns can also affect drug-drug interactions. It is hoped that our work will pave the way for future studies that employ a combination of transcriptomic and complementary techniques to dissect drugdrug interactions.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found at: https://doi.org/10.5061/ dryad.bk3j9kdcn.

AUTHOR CONTRIBUTIONS
QQ and TB conceived the study. SAA and QQ designed, performed, and analyzed the experiments. QQ performed the bioinformatic analyses. TB supervised and managed the project. All authors contributed to the original draft and edited subsequent versions of the manuscript.